# Load libraries, install if needed
pkgs <- c("tidyverse", "readxl", "tidylog", "RCurl", "devtools", "patchwork",
"kableExtra", "viridis", "RColorBrewer", "here", "sdmTMBextra")
if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){
install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
}
invisible(lapply(pkgs, library, character.only = T))
# Packages not on CRAN or dev version
# remotes::install_github("pbs-assess/sdmTMB", dependencies = TRUE)
library(sdmTMB)
# Source code for map plots
# You need: # devtools::install_github("seananderson/ggsidekick") # not on CRAN; library(ggsidekick)
devtools::source_url("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/R/functions/map-plot.R")
options(ggplot2.continuous.colour = "viridis")
# Set path
home <- here::here()Fit tweedie models to biomass density and extract AIC
Intro
Fit Tweedie models to biomass density of cod, flounder, plaice and dab (juveniles and adults) between 1993-2020 in the Baltic Sea using sdmMTB fit different oxygen and temperature-related covariates, compare AIC. Select the “best” covariate for further trend and velocity analysis.
Load packages & source functions
Read data
# Read data
d <- #readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/catch_clean.csv") %>% %>%
readr::read_csv(paste0(home, "/data/clean/catch_clean.csv")) %>%
rename(X = x, Y = y) %>%
pivot_longer(c(cod_adult, cod_juvenile, dab_adult, dab_juvenile, flounder_adult,
flounder_juvenile, plaice_adult, plaice_juvenile),
names_to = "group", values_to = "density") %>%
separate(group, into = c("species", "life_stage"), remove = FALSE) %>%
drop_na(depth, temp, oxy, sal, density)Rows: 12385 Columns: 25
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (4): country, haul_id, ices_rect, month_year
dbl (21): year, lat, lon, quarter, month, sub_div, oxy, temp, sal, depth, x,...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed 2 variables (X, Y)
pivot_longer: reorganized (cod_adult, cod_juvenile, dab_adult, dab_juvenile, plaice_adult, …) into (group, density) [was 12385x25, now 99080x19]
drop_na: removed 9,570 rows (10%), 89,510 rows remaining# Read metabolic parameter estimates and left_join
mi_pars <- #readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/mi_params.csv") %>%
readr::read_csv(paste0(home, "/data/clean/mi_params.csv")) %>%
dplyr::select(n_o2, E_o2, A0_o2, species) # TODO: remove the extra columns for pressure based parametersRows: 4 Columns: 10
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): species
dbl (9): temp, po2, o2, A0_o2, A0_po2, n_po2, E_po2, n_o2, E_o2
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.# TODO: for now we'll use plaice parameters for flounder, see "00_estimate_mi_params.Rmd"
mi_pars# A tibble: 4 × 4
n_o2 E_o2 A0_o2 species
<dbl> <dbl> <dbl> <chr>
1 -0.300 0.294 0.0451 plaice
2 -0.300 0.294 0.0249 cod
3 -0.300 0.294 0.0371 dab
4 -0.300 0.294 0.0120 floundermi_pars <- mi_pars %>%
mutate(A0_o2 = ifelse(species == "flounder",
filter(mi_pars, species == "plaice")$A0_o2,
A0_o2))filter: removed 3 rows (75%), one row remaining
mutate: changed one value (25%) of 'A0_o2' (0 new NA)d <- left_join(d, mi_pars, by = "species")left_join: added 3 columns (n_o2, E_o2, A0_o2)
> rows only in x 0
> rows only in y ( 0)
> matched rows 89,510
> ========
> rows total 89,510# Read size csv to calculate the metabolic index
sizes <- #readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/sizes.csv") %>%
readr::read_csv(paste0(home, "/data/clean/sizes.csv")) %>%
mutate(group = paste(species, name, sep = "_")) %>%
dplyr::select(group, B)Rows: 8 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): species, name
dbl (3): mat_w, max_w, B
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
mutate: new variable 'group' (character) with 8 unique values and 0% NAd <- left_join(d, sizes, by = "group")left_join: added one column (B)
> rows only in x 0
> rows only in y ( 0)
> matched rows 89,510
> ========
> rows total 89,510# Drop dab!
# d %>%
# filter(group == "dab_juvenile" & quarter == 1) %>%
# mutate(pres = ifelse(density > 0, "1", "0")) %>%
# ggplot(aes(X, Y, color = pres)) +
# geom_point(size = 0.3) +
# coord_fixed() +
# facet_wrap(~year)
#
# d %>%
# filter(group == "dab_juvenile" & quarter == 4) %>%
# mutate(pres = ifelse(density > 0, "1", "0")) %>%
# ggplot(aes(X, Y, color = pres)) +
# geom_point(size = 0.3) +
# coord_fixed() +
# facet_wrap(~year)
d <- d %>% filter(!species == "dab")filter: removed 23,124 rows (26%), 66,386 rows remainingCalculate the metabolic index
# Oxygen is ml/L, We want micro mol/L. 1 ml/l = 10^3/22.391 = 44.661 micro mol/l
d$oxy_si <- (d$oxy * (10^3)) / 22.391
# Calculate metabolic indices for a given mass and the temperature and oxygen in data
# Line 123 in https://github.com/fate-spatialindicators/SDM_O2/blob/master/code/mi_functions.R
kb <- 0.000086173324 # Boltzmann's constant
t_ref <- 15 # arbitrary reference temperature
# Calculate the metabolic index
d <- d %>%
mutate(inv_temp = (1/(temp + 273.15) - 1/(t_ref + 273.15)),
phi = A0_o2*(B^n_o2)*oxy_si * exp((E_o2/kb)*inv_temp)) %>%
drop_na(phi)mutate: new variable 'inv_temp' (double) with 12,194 unique values and 0% NA new variable 'phi' (double) with 65,602 unique values and 0% NAdrop_na: no rows removedScale variables
d <- d %>%
group_by(group) %>%
mutate(phi_sc = as.vector(scale(phi)),
temp_sc = as.vector(scale(temp)),
temp_sq = temp_sc^2,
oxy_sc = as.vector(scale(oxy)),
oxy_sq = oxy_sc^2,
sal_sc = as.vector(scale(sal)),
depth_sc = as.vector(scale(depth)),
depth_sq = depth_sc*depth_sc) %>%
mutate(quarter_f = as.factor(quarter),
year_f = as.factor(year)) %>%
ungroup()group_by: one grouping variable (group)mutate (grouped): new variable 'phi_sc' (double) with 65,602 unique values and 0% NA new variable 'temp_sc' (double) with 65,576 unique values and 0% NA new variable 'temp_sq' (double) with 65,576 unique values and 0% NA new variable 'oxy_sc' (double) with 65,568 unique values and 0% NA new variable 'oxy_sq' (double) with 65,568 unique values and 0% NA new variable 'sal_sc' (double) with 65,586 unique values and 0% NA new variable 'depth_sc' (double) with 25,364 unique values and 0% NA new variable 'depth_sq' (double) with 25,364 unique values and 0% NAmutate (grouped): new variable 'quarter_f' (factor) with 2 unique values and 0% NA new variable 'year_f' (factor) with 29 unique values and 0% NAungroup: no grouping variables# Seems like a non-linear relationship with temperature, but a quadratic largely does the job??
ggplot(d, aes(temp_sc, density)) +
geom_smooth(method = "lm", color = "grey50", se = FALSE) +
geom_smooth(method = "gam", formula = y~s(x, k=4), color = "steelblue3", se = FALSE) +
geom_smooth(method = "lm", formula = y ~ x + I(x^2), color = "tomato3", se = FALSE) +
facet_wrap(~group, scales = "free")`geom_smooth()` using formula = 'y ~ x'
ggplot(d, aes(oxy_sc, density)) +
geom_smooth(method = "lm", color = "grey50", se = FALSE) +
geom_smooth(method = "gam", formula = y~s(x, k=4), color = "steelblue3", se = FALSE) +
geom_smooth(method = "lm", formula = y ~ x + I(x^2), color = "tomato3", se = FALSE) +
facet_wrap(~group, scales = "free")`geom_smooth()` using formula = 'y ~ x'
ggplot(d, aes(phi_sc, log(density + 1))) +
geom_smooth(method = "lm", color = "grey50", se = FALSE) +
geom_smooth(method = "gam", formula = y~s(x, k=4), color = "steelblue3", se = FALSE) +
geom_smooth(method = "lm", formula = y ~ x + I(x^2), color = "tomato3", se = FALSE) +
facet_wrap(~group, scales = "free")`geom_smooth()` using formula = 'y ~ x'
# Plot mi by species and life stage
pal <- brewer.pal(name = "Dark2", n = 8)[c(7, 5, 3)]
d %>%
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage)) %>%
ggplot(aes(phi, fill = species, color = species)) +
geom_density(alpha = 0.2) +
geom_vline(xintercept = 0, linetype = 2) +
coord_cartesian(expand = 0) +
facet_wrap(~ life_stage, ncol = 1) +
theme_sleek(base_size = 9) +
theme(legend.position = c(0.85, 0.85)) +
labs(x = "Metabolic index (\u03C6)", y = "Density", fill = "Species", color = "Species") +
scale_color_manual(values = rev(pal), name = "") +
scale_fill_manual(values = rev(pal), name = "")mutate: changed 66,386 values (100%) of 'species' (0 new NA) changed 66,386 values (100%) of 'life_stage' (0 new NA)
ggsave(paste0(home, "/figures/supp/mi_histogram.pdf"), width = 11, height = 11, units = "cm", device = cairo_pdf)Fit models and save AIC
Fit models!
data_list_aic <- list()
for(i in unique(d$group)) {
dd <- d %>% filter(group == i)
mesh <- make_mesh(dd, xy_cols = c("X", "Y"), cutoff = 15) # Need specific meshes because not for all hauls did we record all species
print(
ggplot() +
inlabru::gg(mesh$mesh) +
coord_fixed() +
geom_point(aes(X, Y), data = dd, alpha = 0.2, size = 0.5) +
annotate("text", -Inf, Inf, label = paste("n knots = ", mesh$mesh$n), hjust = -0.1, vjust = 2) +
labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(i, "_", " ")))
)
print(i)
# 0. Linear oxygen and temperature, spatially varying quarter
m0 <- sdmTMB(density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + temp_sc + oxy_sc,
data = dd,
mesh = mesh,
family = tweedie(link = "log"),
spatiotemporal = "IID",
# this setting for spatial field(s) works better for majority of groups compare to single sf
spatial = "off",
spatial_varying = ~0 + quarter_f, # if enabled, keep quarter_f in main formula;
time = "year")
print("m0")
sanity(m0)
dd$m0_res <- residuals(m0)
# 1. Linear oxygen and squared temperature, spatially varying quarter
m1 <- sdmTMB(density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + temp_sc + temp_sq + oxy_sc,
data = dd,
mesh = mesh,
family = tweedie(link = "log"),
spatiotemporal = "IID",
spatial = "off",
spatial_varying = ~0 + quarter_f,
time = "year")
print("m1")
sanity(m1)
dd$m1_res <- residuals(m1)
# 2. interaction between linear oxygen and temperature, drop quadratic
m2 <- sdmTMB(density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + temp_sc * oxy_sc,
data = dd,
mesh = mesh,
family = tweedie(link = "log"),
spatiotemporal = "IID",
spatial = "off",
spatial_varying = ~0 + quarter_f,
time = "year")
print("m2")
sanity(m2)
dd$m2_res <- residuals(m2)
# 3. breakpoint oxygen
m3 <- sdmTMB(density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + temp_sc + temp_sq + breakpt(oxy_sc),
data = dd,
mesh = mesh,
family = tweedie(link = "log"),
spatiotemporal = "IID",
spatial = "off",
spatial_varying = ~0 + quarter_f,
time = "year")
print("m3")
sanity(m3)
dd$m3_res <- residuals(m3)
# 4. logistic oxygen
# m4 <- sdmTMB(density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + temp_sc + temp_sq + logistic(oxy_sc),
# data = dd,
# mesh = mesh,
# family = delta_gamma(link1 = "logit", link2 = "log"),
# spatiotemporal = "off",
# spatial = "on",
# time = "year",
# control = sdmTMBcontrol(newton_loops = 2))
# print("m4")
# sanity(m4)
# dd$m4_res <- residuals(m4)
# 5. linear metabolic index
m5 <- sdmTMB(density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + phi_sc,
data = dd,
mesh = mesh,
family = tweedie(link = "log"),
spatiotemporal = "IID",
spatial = "off",
spatial_varying = ~0 + quarter_f,
time = "year")
print("m5")
sanity(m5)
dd$m5_res <- residuals(m5)
# 6. breakpoint metabolic index!
m6 <- sdmTMB(density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + breakpt(phi_sc),
data = dd,
mesh = mesh,
family = tweedie(link = "log"),
spatiotemporal = "IID",
spatial = "off",
spatial_varying = ~0 + quarter_f,
time = "year")
print("m6")
sanity(m6)
dd$m6_res <- residuals(m6)
# logistic metabolic index!
# m7 <- update(m1, density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + logistic(phi_sc))
# print("m7")
# sanity(m7)
# dd$m7_res <- residuals(m7)
# Plot residuals
p1 <- dd %>%
dplyr::select(m0_res, m1_res, m2_res, m3_res, m5_res, m6_res) %>%
pivot_longer(everything()) %>%
ggplot(aes(sample = value)) +
stat_qq(size = 0.75, shape = 21, fill = NA) +
facet_wrap(~name) +
stat_qq_line() +
labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
theme(aspect.ratio = 1)
print(p1)
# Save AIC
data_list_aic[[i]] <- AIC(m0, m1, m2, m3, m5, m6) %>%
tibble::rownames_to_column("model") %>%
mutate(group = i)
}filter: removed 54,050 rows (81%), 12,336 rows remainingLoading required namespace: INLA[1] "cod_adult"
[1] "m0"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m1"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m2"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m3"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m5"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m6"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably largepivot_longer: reorganized (m0_res, m1_res, m2_res, m3_res, m5_res, …) into (name, value) [was 12336x6, now 74016x2]
mutate: new variable 'group' (character) with one unique value and 0% NAfilter: removed 54,050 rows (81%), 12,336 rows remaining
[1] "cod_juvenile"
[1] "m0"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m1"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m2"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m3"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m5"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m6"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably largepivot_longer: reorganized (m0_res, m1_res, m2_res, m3_res, m5_res, …) into (name, value) [was 12336x6, now 74016x2]
mutate: new variable 'group' (character) with one unique value and 0% NAfilter: removed 54,802 rows (83%), 11,584 rows remaining
[1] "plaice_adult"
[1] "m0"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m1"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m2"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m3"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m5"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m6"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably largepivot_longer: reorganized (m0_res, m1_res, m2_res, m3_res, m5_res, …) into (name, value) [was 11584x6, now 69504x2]
mutate: new variable 'group' (character) with one unique value and 0% NAfilter: removed 54,802 rows (83%), 11,584 rows remaining
[1] "plaice_juvenile"
[1] "m0"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m1"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m2"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m3"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m5"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m6"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably largepivot_longer: reorganized (m0_res, m1_res, m2_res, m3_res, m5_res, …) into (name, value) [was 11584x6, now 69504x2]
mutate: new variable 'group' (character) with one unique value and 0% NAfilter: removed 57,114 rows (86%), 9,272 rows remaining
[1] "flounder_adult"
[1] "m0"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m1"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m2"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m3"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m5"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m6"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably largepivot_longer: reorganized (m0_res, m1_res, m2_res, m3_res, m5_res, …) into (name, value) [was 9272x6, now 55632x2]
mutate: new variable 'group' (character) with one unique value and 0% NAfilter: removed 57,112 rows (86%), 9,274 rows remaining
[1] "flounder_juvenile"
[1] "m0"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m1"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m2"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m3"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m5"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large[1] "m6"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably largepivot_longer: reorganized (m0_res, m1_res, m2_res, m3_res, m5_res, …) into (name, value) [was 9274x6, now 55644x2]
mutate: new variable 'group' (character) with one unique value and 0% NA
# Save aic as data frames
data_aic <- bind_rows(data_list_aic)
write_csv(data_aic, paste0(home, "/output/data_aic_01.csv"))Plot meshes in a single plot using patchwork
g1 <- d %>% filter(group == unique(d$group)[1])filter: removed 54,050 rows (81%), 12,336 rows remainingmesh_g1 <- make_mesh(g1, xy_cols = c("X", "Y"), cutoff = 15)
g2 <- d %>% filter(group == unique(d$group)[5])filter: removed 57,114 rows (86%), 9,272 rows remainingmesh_g2 <- make_mesh(g2, xy_cols = c("X", "Y"), cutoff = 15)
g3 <- d %>% filter(group == unique(d$group)[3])filter: removed 54,802 rows (83%), 11,584 rows remainingmesh_g3 <- make_mesh(g3, xy_cols = c("X", "Y"), cutoff = 15)
g4 <- d %>% filter(group == unique(d$group)[2])filter: removed 54,050 rows (81%), 12,336 rows remainingmesh_g4 <- make_mesh(g4, xy_cols = c("X", "Y"), cutoff = 15)
g5 <- d %>% filter(group == unique(d$group)[6])filter: removed 57,112 rows (86%), 9,274 rows remainingmesh_g5 <- make_mesh(g5, xy_cols = c("X", "Y"), cutoff = 15)
g6 <- d %>% filter(group == unique(d$group)[4])filter: removed 54,802 rows (83%), 11,584 rows remainingmesh_g6 <- make_mesh(g6, xy_cols = c("X", "Y"), cutoff = 15)
p1 <- ggplot() +
inlabru::gg(mesh_g1$mesh) +
coord_fixed() +
geom_point(aes(X, Y), data = g1, alpha = 0.2, size = 0.5) +
annotate("text", -Inf, Inf, label = paste("n knots = ", mesh_g1$mesh$n), hjust = -0.1, vjust = 2, size = 2.5) +
labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(g1$group[1], "_", " "))) +
theme_sleek(base_size = 8) +
theme(axis.title.x = element_blank(),
plot.margin = unit(c(0,0,0,0), "cm"))
p2 <- ggplot() +
inlabru::gg(mesh_g2$mesh) +
coord_fixed() +
geom_point(aes(X, Y), data = g2, alpha = 0.2, size = 0.5) +
annotate("text", -Inf, Inf, label = paste("n knots = ", mesh_g2$mesh$n), hjust = -0.1, vjust = 2, size = 2.5) +
labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(g2$group[], "_", " "))) +
theme_sleek(base_size = 8) +
theme(axis.title.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.title.y = element_blank(),
plot.margin = unit(c(0,0,0,0), "cm"))
p3 <- ggplot() +
inlabru::gg(mesh_g3$mesh) +
coord_fixed() +
geom_point(aes(X, Y), data = g3, alpha = 0.2, size = 0.5) +
annotate("text", -Inf, Inf, label = paste("n knots = ", mesh_g3$mesh$n), hjust = -0.1, vjust = 2, size = 2.5) +
labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(g3$group[], "_", " "))) +
theme_sleek(base_size = 8) +
theme(axis.title.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.title.y = element_blank(),
plot.margin = unit(c(0,0,0,0), "cm"))
p4 <- ggplot() +
inlabru::gg(mesh_g4$mesh) +
coord_fixed() +
geom_point(aes(X, Y), data = g4, alpha = 0.2, size = 0.5) +
annotate("text", -Inf, Inf, label = paste("n knots = ", mesh_g4$mesh$n), hjust = -0.1, vjust = 2, size = 2.5) +
labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(g4$group[], "_", " "))) +
theme_sleek(base_size = 8) +
theme(axis.title.x = element_blank(),
plot.margin = unit(c(0,0,0,0), "cm"))
p5 <- ggplot() +
inlabru::gg(mesh_g5$mesh) +
coord_fixed() +
geom_point(aes(X, Y), data = g5, alpha = 0.2, size = 0.5) +
annotate("text", -Inf, Inf, label = paste("n knots = ", mesh_g5$mesh$n), hjust = -0.1, vjust = 2, size = 2.5) +
labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(g5$group[], "_", " "))) +
theme_sleek(base_size = 8) +
theme(axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
plot.margin = unit(c(0,0,0,0), "cm"))
p6 <- ggplot() +
inlabru::gg(mesh_g6$mesh) +
coord_fixed() +
geom_point(aes(X, Y), data = g6, alpha = 0.2, size = 0.5) +
annotate("text", -Inf, Inf, label = paste("n knots = ", mesh_g6$mesh$n), hjust = -0.1, vjust = 2, size = 2.5) +
labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(g6$group[], "_", " "))) +
theme_sleek(base_size = 8) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
plot.margin = unit(c(0,0,0,0), "cm"))
p1 + p2 + p3 + p4 + p5 + p6
ggsave(paste0(home, "/figures/supp/mesh.pdf"), width = 17, height = 11, units = "cm", device = cairo_pdf)